function [psicf20, rhojcf20, sjcf20, sHjcf20,sNHjcf20,  s0jcf20, s0Hjcf20,s0NHjcf20, dealer_idscf20,qcf20,etaDjcf20,etaEjcf20,uBcf20,uEcf20,piBcf20,piEcf20,YBcf20,YEcf20,Ecostcf20,uE2cf20, W20, Wxij20, Wvdj20, WxijH20, W20_shrink, W20wH, Wxij20wH, Wvdj20wH, WxijH20wH, W20_shrinkwH, psicf20_g,rhojcf20_g,sjcf20_g,sHjcf20_g, sNHjcf20_g, s0jcf20_g, s0Hjcf20_g, s0NHjcf20_g, uBcf20_g,uEcf20_g,YBcf20_g,YEcf20_g,Ecostcf20_g, piBcf20_g,piEcf20_g, uE2cf20_g, psicf20_opt, rhojcf20_opt, psicf20_opt_g, rhojcf20_opt_g, W20_g, Wxij20_g, Wvdj20_g, WxijH20_g, W20_shrink_g, W20_gwH, Wxij20_gwH, Wvdj20_gwH, WxijH20_gwH, W20_shrink_gwH] = counterfac_per_side_fun(CF, side, access, Y, new_period, new_period_side, vdjSQ, kappa, flexible)
  
global T Imax Dmax homedealer eta loyalty

%Set option for maximization
options = optimoptions('fmincon','FiniteDifferenceType','central','Algorithm', 'interior-point', 'OptimalityTolerance',1e-40, 'StepTolerance', 1e-40, 'FunctionTolerance', 1e-40,   'ConstraintTolerance', 1e-40, 'MaxFunctionEvaluations', 10^6, 'MaxIterations', 10^6);  
lb =[0 0 0 0 0 0 0 0 0];

%Initialize (flexible quotes at estimated cost) 
psicf20     = NaN* ones(T,Imax);
rhojcf20    = NaN* ones(T,Imax);
psicf20_opt = NaN* ones(T,Imax);
rhojcf20_opt= NaN* ones(T,Imax); 
sjcf20      = NaN* ones(T,Imax); 
sHjcf20     = NaN* ones(T,Imax); 
sNHjcf20    = NaN* ones(T,Imax, Imax); 
s0Hjcf20     = NaN* ones(T,Imax); 
s0NHjcf20    = NaN* ones(T,Imax, Imax); 
s0jcf20     = NaN* ones(T,Imax); 
dealer_idscf20=NaN* ones(T,Imax);
qcf20       = NaN* ones(T,Imax);
etaDjcf20   = NaN* ones(T,Imax); 
etaEjcf20   = NaN* ones(T,Imax); 
uBcf20      = NaN* ones(T,Imax); 
uEcf20      = NaN* ones(T,Imax); 
uE2cf20     = NaN* ones(T,Imax); 
piBcf20     = NaN* ones(T,Imax); 
piEcf20     = NaN* ones(T,Imax); 
YBcf20      = NaN* ones(T,Imax); 
YEcf20      = NaN* ones(T,Imax); 
Ecostcf20   = NaN* ones(T,Imax);
W20         = NaN* ones(T,1);
W20_shrink  = NaN* ones(T,1);
Wxij20      = NaN* ones(T,1);
Wvdj20      = NaN* ones(T,1);
WxijH20     = NaN* ones(T,1);
W20wH         = NaN* ones(T,1);
W20_shrinkwH  = NaN* ones(T,1);
Wxij20wH      = NaN* ones(T,1);
Wvdj20wH      = NaN* ones(T,1);
WxijH20wH     = NaN* ones(T,1);


%Initialize (flexible quotes at estimated cost) for I
psicf20_g     = NaN* ones(T,2*Imax);
rhojcf20_g    = NaN* ones(T,2*Imax); 
psicf20_opt_g  = NaN* ones(T,2*Imax);
rhojcf20_opt_g = NaN* ones(T,2*Imax); 
sjcf20_g      = NaN* ones(T,2*Imax); 
sHjcf20_g     = NaN* ones(T,2*Imax); 
sNHjcf20_g    = NaN* ones(T,2*Imax, 2*Imax); 
s0Hjcf20_g     = NaN* ones(T,2*Imax); 
s0NHjcf20_g    = NaN* ones(T,2*Imax, 2*Imax); 
s0jcf20_g     = NaN* ones(T,2*Imax);  
uBcf20_g      = NaN* ones(T,2*Imax); 
uEcf20_g      = NaN* ones(T,2*Imax); 
uE2cf20_g     = NaN* ones(T,2*Imax); 
YBcf20_g      = NaN* ones(T,2*Imax); 
YEcf20_g      = NaN* ones(T,2*Imax); 
Ecostcf20_g   = NaN* ones(T,2*Imax);
piBcf20_g     = NaN* ones(T,2*Imax);
piEcf20_g     = NaN* ones(T,2*Imax);
W20_g         = NaN* ones(T,2);
Wxij20_g      = NaN* ones(T,2);
Wvdj20_g      = NaN* ones(T,2);
WxijH20_g     = NaN* ones(T,2);
W20_shrink_g  = NaN* ones(T,2);
W20_gwH         = NaN* ones(T,2);
Wxij20_gwH      = NaN* ones(T,2);
Wvdj20_gwH      = NaN* ones(T,2);
WxijH20_gwH     = NaN* ones(T,2);
W20_shrink_gwH  = NaN* ones(T,2);

for k=1:T
   
    date=new_period(k,1);
    new_period_now = new_period_side(new_period_side(:,4)==side,1:3);
    
    t=find(new_period_now(:,1)==date); %line in new-period of the side of the trade

    %If this data was a trading date on that size
    if ~isempty(t) 

    idxt = [new_period_now(t,2):new_period_now(t,3)]'; %lines in Y
 
    %label
    Yt = Y(idxt,:);
       
    if side==1
        
        paramst_I = Yt(Yt(:,1)==1,28:30); %IB
        paramst_R = Yt(Yt(:,1)==2,28:30); %RB

        %both groups were active
        if sum(Yt(:,1)==1)>0  && sum(Yt(:,1)==2)>0
            paramst  = [paramst_I(1,:);paramst_R(1,:)];
        %only retailers active (unse average of estimated parameters)    
        elseif  sum(Yt(:,1)==1)==0  && sum(Yt(:,1)==2)>0
            paramst  = [nanmean( Y(Y(:,1)==1,28:30));paramst_R(1,:)];
        %only institutionals active (unse average of estimated parameters)    
        elseif sum(Yt(:,1)==1)>0  && sum(Yt(:,1)==2)==0
            paramst  = [paramst_I(1,:);nanmean( Y(Y(:,1)==2,28:30))];
        end 

    elseif side==-1
            
        paramst_I = Yt(Yt(:,1)==3,28:30); %IS
        paramst_R = Yt(Yt(:,1)==4,28:30); %RS

        if sum(Yt(:,1)==3)>0  && sum(Yt(:,1)==4)>0
            paramst   = [paramst_I(1,:);paramst_R(1,:)];
        elseif sum(Yt(:,1)==3)==0  && sum(Yt(:,1)==4)>0
            paramst  = [nanmean(Y(Y(:,1)==3,28:30));paramst_R(1,:)];
        elseif sum(Yt(:,1)==3)>0  && sum(Yt(:,1)==4)==0
            paramst  = [paramst_I(1,:);nanmean(Y(Y(:,1)==4,28:30))];
        end 
    end 
     

    %If cost is 0 in counterfact
    if CF==2|| CF==4 || CF==11 || CF==15 || CF==6 || CF==10 || CF==14 , paramst(:,3)=[0;0];end %eliminate all costs
    if CF==5 || CF==9 || CF==13 , paramst(:,3)=paramst(:,3)+1.1; end %eliminate the fee
   
    if (CF==6 || CF==10 || CF==14 ) || loyalty==0 %set home dealer benefit to zero
        Yt(:,25) =Yt(:,24);
    end 

    if homedealer==0 || (homedealer==1 && CF~=7) 
        params_combinedt = (1-kappa).*paramst(1,:)+ kappa.*paramst(2,:);
    elseif homedealer==1 && CF==7 %only institutionals
        params_combinedt = paramst(1,:);
    end 
    
    %If parameters were estimated for all groups and dealer's values are estimated
    if ~isnan(sum(params_combinedt)) && nansum(vdjSQ(t,:))>0
 
    if homedealer==0
       
       %index for dealers
        l=1;
        for j=1:length(Dmax)
            if sum(Yt(:,17)==Dmax(j))>0
            idxj(l)=find(Yt(:,17)==Dmax(j), 1, 'first');
            l=l+1;
            end
        end

        q0      = Yt(idxj, 19);     
        vdjSQt  = vdjSQ(t,:);
        
        %For CF's where q=v
        if CF>=8 && CF<=15 ||  (CF==17 || CF==18) , q0=vdjSQt'; end 
        
     if eta==0
         
        if flexible==1 %find quotes if quotes are flexible 
            [qcf20t,fval(t,:),flag(t,:),output] = fmincon(@(q) model_prediction_avg(q, Yt(:,2:end), params_combinedt, side,0, CF,access, vdjSQt,idxj,t), q0, [],[],[],[],lb,[],[],options);        
        else %use observed quotes if quotes are fixed
            qcf20t = q0;
        end 

        %All output with representative investor 
        [OFcf20t(t,:), psicf20t, rhojcf20t, sjcf20t, s0jcf20t,etaEj20t,etaDj20t, vdj20t, piDprimej20t,  dealer_idscf20t, uBcf20t, uEcf20t,  piBcf20t, piEcf20t, YBcf20t, YEcf20t, Ecostcf20t, uE2cf20t, psicf20_optt, rhojcf20_optt, W20t, Wxij20t, Wvdj20t, W20_schrinkt] = model_prediction_avg(qcf20t, Yt(:,2:end), params_combinedt, side, 0, CF, access, vdjSQt,idxj, t);

        %All output for each investor group separately
        [a0(t,:), psicf20t_I, rhojcf20t_I, sjcf20t_I, s0jcf20t_I,etaEj20t_I,etaDj20t_I, vdj20t_I, piDprimej20t_I,  dealer_idscf20t_I, uBcf20t_I, uEcf20t_I,  piBcf20t_I, piEcf20t_I, YBcf20t_I, YEcf20t_I, Ecostcf20t_I, uE2cf20t_I, psicf20_optt_I, rhojcf20_optt_I, W20t_I, Wxij20t_I, Wvdj20t_I, W20_schrinkt_I] = model_prediction_avg(qcf20t, Yt(:,2:end), paramst(1,:), side, 0, CF, access, vdjSQt,idxj, t);
        [a1(t,:), psicf20t_R, rhojcf20t_R, sjcf20t_R, s0jcf20t_R,etaEj20t_R,etaDj20t_R, vdj20t_R, piDprimej20t_R,  dealer_idscf20t_R, uBcf20t_R, uEcf20t_R,  piBcf20t_R, piEcf20t_R, YBcf20t_R, YEcf20t_R, Ecostcf20t_R, uE2cf20t_R, psicf20_optt_R, rhojcf20_optt_R, W20t_R, Wxij20t_R, Wvdj20t_R, W20_schrinkt_R] = model_prediction_avg(qcf20t, Yt(:,2:end), paramst(2,:), side, 0, CF, access, vdjSQt,idxj, t);
   
     elseif eta>0 
         
        if flexible==1 %find quotes if quotes are flexible 
            [qcf20t,fval(t,:),flag(t,:),output] = fmincon(@(q) model_prediction_avg_eta(q, Yt(:,2:end), params_combinedt, side,0,access, vdjSQt,idxj,CF), q0, [],[],[],[],lb,[],[],options);        
        else %use observed quotes if quotes are fixed
            qcf20t = q0;
        end 

        %All output with representative investor 
        [OFcf20t(t,:), psicf20t, rhojcf20t, sjcf20t, s0jcf20t,etaEj20t,etaDj20t, vdj20t, piDprimej20t,  dealer_idscf20t, uBcf20t, uEcf20t,  piBcf20t, piEcf20t, YBcf20t, YEcf20t, Ecostcf20t, uE2cf20t, psicf20_optt, rhojcf20_optt, W20t, Wxij20t, Wvdj20t, W20_schrinkt] = model_prediction_avg_eta(qcf20t, Yt(:,2:end), params_combinedt, side,0, access, vdjSQt,idxj, CF);

        %All output for each investor group separately
        [a0(t,:), psicf20t_I, rhojcf20t_I, sjcf20t_I, s0jcf20t_I,etaEj20t_I,etaDj20t_I, vdj20t_I, piDprimej20t_I,  dealer_idscf20t_I, uBcf20t_I, uEcf20t_I,  piBcf20t_I, piEcf20t_I, YBcf20t_I, YEcf20t_I, Ecostcf20t_I, uE2cf20t_I, psicf20_optt_I, rhojcf20_optt_I, W20t_I, Wxij20t_I, Wvdj20t_I, W20_schrinkt_I] = model_prediction_avg_eta(qcf20t, Yt(:,2:end), paramst(1,:), side, 0, access, vdjSQt,idxj, CF);
        [a1(t,:), psicf20t_R, rhojcf20t_R, sjcf20t_R, s0jcf20t_R,etaEj20t_R,etaDj20t_R, vdj20t_R, piDprimej20t_R,  dealer_idscf20t_R, uBcf20t_R, uEcf20t_R,  piBcf20t_R, piEcf20t_R, YBcf20t_R, YEcf20t_R, Ecostcf20t_R, uE2cf20t_R, psicf20_optt_R, rhojcf20_optt_R, W20t_R, Wxij20t_R, Wvdj20t_R, W20_schrinkt_R] = model_prediction_avg_eta(qcf20t, Yt(:,2:end), paramst(2,:), side, 0, access, vdjSQt,idxj, CF);
       
     end 
        
   elseif homedealer==1 && eta==0
                     
        %index for dealers
        l=1;
        for j=1:length(Dmax)
            if sum(Yt(:,17)==Dmax(j))>0
            idxj(l)=find(Yt(:,17)==Dmax(j), 1, 'first');
            l=l+1;
            end
        end
    
        q0      = Yt(idxj, 19);
        vdjSQt  = vdjSQ(t,:);
        
        %For CF's where q=v
        if CF>=8 && CF<=15  || (CF==17 || CF==18), q0=vdjSQt'; end 
        
        if flexible==1 %find quotes if quotes are flexible 
            [qcf20t,fval(t,:),flag(t,:),output] = fmincon(@(q) model_prediction_avg_H(q, Yt(:,2:end), params_combinedt, side,0, CF,access, vdjSQ(t,:), idxj), q0, [],[],[],[],lb,[],[],options);        
        else %use observed quotes if quotes are fixed
            qcf20t = q0;
        end 

        %All output with representative investor 
        [OFcf20t(t,:), psicf20t, rhojcf20t, sjcf20t, sHjcf20t, sNHjcf20t, s0jcf20t, s0Hjcf20t, s0NHjcf20t, etaEj20t,etaDj20t, vdj20t, piDprimej20t,  dealer_idscf20t, uBcf20t, uEcf20t,  piBcf20t, piEcf20t, YBcf20t, YEcf20t, Ecostcf20t, uE2cf20t, psicf20_optt, rhojcf20_optt, ...
            W20t, Wxij20t, Wvdj20t, WxijH20t, W20_schrinkt, W20twH, Wxij20twH, Wvdj20twH, WxijH20twH, W20_schrinktwH ] = model_prediction_avg_H(qcf20t, Yt(:,2:end), params_combinedt, side, 0, CF, access, vdjSQ(t,:), idxj);

        %All output for each investor group separately
        [a0(t,:), psicf20t_I, rhojcf20t_I, sjcf20t_I, sHjcf20t_I, sNHjcf20t_I, s0jcf20t_I, s0Hjcf20t_I, s0NHjcf20t_I, etaEj20t_I,etaDj20t_I, vdj20t_I, piDprimej20t_I,  dealer_idscf20t_I, uBcf20t_I, uEcf20t_I,  piBcf20t_I, piEcf20t_I, YBcf20t_I, YEcf20t_I, Ecostcf20t_I, uE2cf20t_I, psicf20_optt_I, rhojcf20_optt_I, ...
            W20t_I, Wxij20t_I, Wvdj20t_I, WxijH20t_I, W20_schrinkt_I, W20t_IwH, Wxij20t_IwH, Wvdj20t_IwH, WxijH20t_IwH, W20_schrinkt_IwH] = model_prediction_avg_H(qcf20t, Yt(:,2:end), paramst(1,:), side, 0, CF, access, vdjSQ(t,:), idxj);
        
        %If retail investors earn no home dealer benefit 
        YtR = Yt;
        [a1(t,:), psicf20t_R, rhojcf20t_R, sjcf20t_R, sHjcf20t_R, sNHjcf20t_R, s0jcf20t_R, s0Hjcf20t_R, s0NHjcf20t_R, etaEj20t_R,etaDj20t_R, vdj20t_R, piDprimej20t_R,  dealer_idscf20t_R, uBcf20t_R, uEcf20t_R,  piBcf20t_R, piEcf20t_R, YBcf20t_R, YEcf20t_R, Ecostcf20t_R, uE2cf20t_R, psicf20_optt_R, rhojcf20_optt_R, ...
            W20t_R, Wxij20t_R, Wvdj20t_R, WxijH20t_R, W20_schrinkt_R, W20t_RwH, Wxij20t_RwH, Wvdj20t_RwH, WxijH20t_RwH, W20_schrinkt_RwH] = model_prediction_avg_H(qcf20t, YtR(:,2:end), paramst(2,:), side, 0, CF, access, vdjSQ(t,:), idxj);
   end


   %Store and stack  
    idxd              = find(ismember(Dmax,dealer_idscf20t));      
    qcf20t_full       = NaN.*ones(Imax,1);
    if CF>=8 && CF<=15 || (CF==17 || CF==18)
        qcf20t_full(:) = qcf20t; 
    else 
        qcf20t_full(idxd) = qcf20t;
    end 
    qcf20(k,:)          = qcf20t_full'; 
    etaDjcf20(k,:)      = etaDj20t';
    etaEjcf20(k,:)      = etaEj20t';
    psicf20(k,:)        = psicf20t'; 
    rhojcf20_opt(k,:)   = rhojcf20_optt';
    psicf20_opt(k,:)    = psicf20_optt'; 
    rhojcf20(k,:)       = rhojcf20t';
    sjcf20(k,:)         = sjcf20t';
    
    if homedealer==1
    	sHjcf20(k,:)        = sHjcf20t';
        sNHjcf20(k,:,:)     = sNHjcf20t';
        s0Hjcf20(k,:)       = s0Hjcf20t';
        s0NHjcf20(k,:,:)    = s0NHjcf20t';

    end
    
    s0jcf20(k,:)        = s0jcf20t';
    dealer_idscf20(k,:) = dealer_idscf20t';
    uBcf20(k,:)         = uBcf20t';
    uEcf20(k,:)         = uEcf20t';
    uE2cf20(k,:)        = uE2cf20t';
    piBcf20(k,:)        = piBcf20t';
    piEcf20(k,:)        = piEcf20t';
    YBcf20(k,:)         = YBcf20t';
    YEcf20(k,:)         = YEcf20t';
    Ecostcf20(k,:)      = Ecostcf20t';
    W20(k,:)            = W20t';
    W20_shrink(k,:)     = W20_schrinkt';
    Wxij20(k,:)         = Wxij20t';
    Wvdj20(k,:)         = Wvdj20t';
    
    if homedealer==1
        WxijH20(k,:)      = WxijH20t';
        W20wH(k,:)        = W20twH';
        W20_shrinkwH(k,:) = W20_schrinktwH';
        Wxij20wH(k,:)     = Wxij20twH';
        Wvdj20wH(k,:)     = Wvdj20twH';
        WxijH20wH(k,:)    = WxijH20twH';
    end 
    
    psicf20_g(k,1:Imax)    = psicf20t_I';
    rhojcf20_g(k,1:Imax)   = rhojcf20t_I'; 
    psicf20_opt_g(k,1:Imax) = psicf20_optt_I';
    rhojcf20_opt_g(k,1:Imax)= rhojcf20_optt_I';
    sjcf20_g(k,1:Imax)     = sjcf20t_I'; 
    
    if homedealer==1
        sHjcf20_g(k,1:Imax)              = sHjcf20t_I';  
        sNHjcf20_g(k,1:Imax, 1:Imax)     = sNHjcf20t_I';  
        s0Hjcf20_g(k,1:Imax)             = s0Hjcf20t_I';  
        s0NHjcf20_g(k,1:Imax, 1:Imax)    = s0NHjcf20t_I';  

    end 
    s0jcf20_g(k,1:Imax)    = s0jcf20t_I';  
    uBcf20_g(k,1:Imax)     = uBcf20t_I'; 
    uEcf20_g(k,1:Imax)     = uEcf20t_I'; 
    uE2cf20_g(k,1:Imax)    = uE2cf20t_I'; 
    YBcf20_g(k,1:Imax)     = YBcf20t_I'; 
    YEcf20_g(k,1:Imax)     = YEcf20t_I'; 
    Ecostcf20_g(k,1:Imax)  = Ecostcf20t_I';
    piBcf20_g(k,1:Imax)    = piBcf20t_I'; 
    piEcf20_g(k,1:Imax)    = piEcf20t_I'; 
    
    W20_g(k,1)            = W20t_I';
    W20_shrink_g(k,1)     = W20_schrinkt_I';
    Wxij20_g(k,1)         = Wxij20t_I';
    Wvdj20_g(k,1)         = Wvdj20t_I';
    
    if homedealer==1
       WxijH20_g(k,1)          = WxijH20t_I';
       W20_gwH(k,1)            = W20t_IwH';
       W20_shrink_gwH(k,1)     = W20_schrinkt_IwH';
       Wxij20_gwH(k,1)         = Wxij20t_IwH';
       Wvdj20_gwH(k,1)         = Wvdj20t_IwH';
       WxijH20_gwH(k,1)        = WxijH20t_IwH';
    end 
    
    
    psicf20_g(k,(Imax+1):end)       = psicf20t_R';
    rhojcf20_g(k,(Imax+1):end)      = rhojcf20t_R'; 
    psicf20_opt_g(k,(Imax+1):end)   = psicf20_optt_R';
    rhojcf20_opt_g(k,(Imax+1):end)  = rhojcf20_optt_R'; 
    sjcf20_g(k,(Imax+1):end)        = sjcf20t_R'; 
    
    if homedealer==1
        sHjcf20_g(k,(Imax+1):end)                = sHjcf20t_R';  
        sNHjcf20_g(k,(Imax+1):end,(Imax+1):end)  = sNHjcf20t_R';
        
        s0Hjcf20_g(k,(Imax+1):end)                = s0Hjcf20t_R';  
        s0NHjcf20_g(k,(Imax+1):end,(Imax+1):end)  = s0NHjcf20t_R';

    end 
    s0jcf20_g(k,(Imax+1):end)  = s0jcf20t_R';  
    uBcf20_g(k,(Imax+1):end)   = uBcf20t_R'; 
    uEcf20_g(k,(Imax+1):end)   = uEcf20t_R'; 
    uE2cf20_g(k,(Imax+1):end)  = uE2cf20t_R'; 
    YBcf20_g(k,(Imax+1):end)   = YBcf20t_R'; 
    YEcf20_g(k,(Imax+1):end)   = YEcf20t_R'; 
    Ecostcf20_g(k,(Imax+1):end)  = Ecostcf20t_R';
    piBcf20_g(k,(Imax+1):end)    = piBcf20t_R'; 
    piEcf20_g(k,(Imax+1):end)    = piEcf20t_R'; 
    
    W20_g(k,2)            = W20t_R';
    W20_shrink_g(k,2)     = W20_schrinkt_R';
    Wxij20_g(k,2)         = Wxij20t_R';
    Wvdj20_g(k,2)         = Wvdj20t_R';
    
    if homedealer==1
       WxijH20_g(k,2)          = WxijH20t_R';
       W20_gwH(k,2)            = W20t_RwH';
       W20_shrink_gwH(k,2)     = W20_schrinkt_RwH';
       Wxij20_gwH(k,2)         = Wxij20t_RwH';
       Wvdj20_gwH(k,2)         = Wvdj20t_RwH';
       WxijH20_gwH(k,2)        = WxijH20t_RwH';
    end 

    clear idxd idxj
    end 
    end 
 
end 